# sandwich estimator function

summarize.r = function(reg.lm, ...)
{
  z = summary.lm(reg.lm)
  W.VCOV = vcovHC(reg.lm,type="HC0")
  sig.sq = ((summary(reg.lm)$sigma)^2)
  coef = z$coefficients[,1]
  
  White.SE <- sqrt(diag(W.VCOV))
  t.robust <- coef/White.SE ## Element-by-element division
  df <- reg.lm$df.residual
  p.val.robust <- 2*(1-pt(abs(t.robust),df))
  z$cov.unscaled = W.VCOV/sig.sq
  z$coefficients[,2] = White.SE
  z$coefficients[,3] = t.robust
  z$coefficients[,4] = p.val.robust

  return(z)
}
